Negative velocity fluctuations of pulled reaction fronts 
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The position of a reaction front, propagating into an unstable state, fluctuates because of the 
shot noise. What is the probability that the fluctuating front moves considerably slower than its 
deterministic counterpart? Can the noise arrest the front motion for some time, or even make it 
move in the wrong direction? We present a WKB theory that assumes many particles in the front 
region and answers these questions for the microscopic model A £2 2A and random walk. 
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The Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) 
equation 



d t q = q-q 2 + d 2 x q , 



(1) 



describes invasion of an unstable phase, q(x — > oo, t) = 0, 
by a stable phase, q(x — > — oo, t) = 1. It is one of the 
most fundamental models in mathematical genetics and 
population biology 111, |2[ , but similar equations appear in 
chemical kinetics U extreme value statistics ji], disor- 
dered systems [B| and even particle physics [||. 

The invasion fronts are described by traveling-front 
solutions (TFSs) of Eq. (5): q(x,t) = Q , C (Q, where 
£ = x — ct. Qo,c(0 solves the equation 



Oo, c + cQo,c + Qo,c-Q2, c = : 



(2) 



where the prime denotes the ^-derivative. It is well 
known [3] that, for a steep enough initial condition, the 
invasion front of Eq. |T]) converges at long times to the 
limiting TFS, Qo,a(0> of Ec l- © with the velocity c = 2. 
This special velocity is determined by the dynamics of the 
leading edge of the front, where Eq. ([I]) can be linearized 
around q = 0. One can say that the nonlinear front, de- 
scribed by Eq. (fTJ) , is "pulled" by its leading edge, hence 
the term "pulled fronts" Q , of which the FKPP equation 
((T|) is the celebrated example. 

Being a mean-field equation, Eq. (p} ignores discrete- 
ness of particles and the resulting shot noise. Their im- 
pact on the front propagation is dramatic, and it has 
attracted a great deal of interest I« jl6 f. The noise leads 
to fluctuations of the front shape [8|, [l2| . The particle 
discreteness and noise also cause deviations of the front 
position with time that include a systematic part - the 
front velocity shift - and a fluctuating part. If JV ^> 1 
is the number of particles in the front region, the front 
velocity shift scales as In -2 N @-ll|, whereas the front 
diffusion coefficient scales as In -3 N 

[HE EH. These 
properties of fluctuating pulled fronts are markedly dif- 
ferent from those of fluctuating fronts propagating into 
metastable states, where the front velocity shift and the 



front diffusion coefficient scale as inverse powers of N and 
are therefore much smaller 



1J, 113,118 1 . 



The front diffusion coefficient probes typical, relatively 
small fluctuations of the front position. Here we ask the 
following question that has not been addressed before. 
What is the probability V(c) that a fluctuating front 
moves, during a long time interval r, with average ve- 
locity c that is considerably smaller than the mean-field 
value Co = 2? This includes the extreme case of c = 0, 
when the front is standing on average, and even c < 0, 
when it moves "in the wrong direction." We will also 
obtain new results in the regime of 2 — c <C 1. Impor- 
tantly, these results can be extended to a broad class of 
pulled fronts. Still, we will focus on a specific microscopic 
model which mean- field limit is Eq. ([T]). A WKB (after 
Wentzcl, Kramers and Brillouin) formalism transforms 
the governing master equation into an effective Hamilto- 
nian mechanics, continuous in space and time. Then the 
problem can be reduced to finding a TFS, with given c, 
of the Hamilton's equations. This TFS turns out to be 
instanton-like, as it connects two "equilibria." The action 
along the instanton yields, up to a pre-factor, Vic) which 
is strongly non-Gaussian, exponentially small in N and 
rapidly falls as c goes down. 

The microscopic model that we adopt here, see e.g. 



Ref. [ljj ll2|, involves particles on a one-dimensional 
lattice with lattice constant h. The particles undergo 
stochastic on-site reactions A ~ > 2A and 2A — > A with 
rate constants A and 2X/K, respectively, and unbiased 
random walk between neighboring sites with rate con- 
stant Dq. For K ^> 1 and Dq ^> A, the mean- field theory 
for this problem is Eq. ([T]), see e.g. Ref. [19|. Here q is 
(the continuous limit of) the on-site number of particles 
rii rescaled by K, time is rescaled by 1/ A, a; is rescaled by 
Id = (-D/A) 1 / 2 , and D = Doh 2 is the diffusion coefficient. 
The front velocity is measured in units of (AD) 1 / 2 . 

The stochastic problem is formulated in terms of the 
master equation for the multivariate probability distribu- 
tion 7r(n; t) fiol. f]~8l - j2oj | which describes how the probabil- 
ity of having rii particles on site i at time t changes via 
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each of the elementary processes. We apply the WKB 
ansatz 7r(n;t) = exp[— KS(q\ t)], with q = n/K, to 
the master equation, see e.g. Ref. (l9j . As K S> 1, 
we can treat S as a smooth functional that can be 
Taylor expanded in q. Furthermore, in the fast dif- 
fusion regime, Dq 3> A, we can use a continuous de- 
scription in x. In the leading order in 1/K, the master 
equation reduces to a functional Hamilton- Jacobi equa- 
tion dtS + H(q, dqS) = 0. The Hamiltonian functional 
H = h^ 1 J dx w, where 

w = q(e p - 1) + q 2 (e-P - 1) - d x qd x p + q(d x p) 2 . (3) 

The Hamilton's equations for the "coordinate" q{x,t) 
(the particle number density) and conjugate momentum 
p{x, t) are 

d t q = qe p - q 2 e-'P + d 2 x q - 2d x (qd x p) , (4) 
dtp = -(e p - 1) - 2q(e-P - 1) - d 2 xP - {d x pf . (5) 

We assume that the initial particle density is zero to 
the right of some finite point in space. At x — > — oo 
there is a stationary distribution of the particle density, 
peaked at q = 1. Therefore, we demand q(—oo,t) = 1 
and p{— oo, f) = which corresponds to the fixed point 
(q = l,p = 0) of the on-site Hamiltonian Ho(q,p) = 
q(e p — 1) + q 2 {e~ p — 1). We also demand q(oo,t) = 0. 
The momentum p can be unbounded at x = oo, but must 
be bounded at finite x. 

The calculations greatly simplify in the new variables 
Q = qe-P and P = e p - 1 [HI, SI]. The generating func- 
tion of this canonical transformation is h^ 1 J dxF{q, Q), 
where 

F(q,Q) = qln(q/Q)-q + Q. (6) 

The new Hamiltonian density becomes 

W = (Q-Q 2 )(P + P 2 )-d x Qd x P, (7) 

whereas the Hamilton's equations are 

d t Q = Q-Q 2 + 2QP-2Q 2 P + dlQ, (8) 
d t P = -P - P 2 + 2QP + 2QP 2 - d 2 x P . (9) 

The mean-field FKPP equation ([1]) corresponds to the 
motion on the invariant manifold P(x,t) = 0. 

The boundary conditions (BCs) at x = ±oo are de- 
termined by the zero-energy fixed points (Q, P) of the 
on-site Hamiltonian H = (Q - Q 2 )(P + P 2 ). There are 
four such points: (1,0), (0,0), (0,-1), and (1,-1). In 
view of the BCs for q and p at x — — oo, we demand 
Q(-oo,t) = 1 and P(-oo,<) = 0. Now, as q = Q(P+1), 
the BC q(oo,t) = can be satisfied at any of the three 
fixed points: Q(oo,i) = 0, P(oo,£) = —1 (we will call 
these BCs "BCs A"), Q(oo,t) = 1, P(oo,f) = -1 (BCs 
B), or trivial BCs Q(oo,i) = P(oo,t) — 0. It turns out 



that each of the BCs A and B is possible for fluctuat- 
ing fronts, whereas the trivial BCs are only possible for 
mean-field fronts, where P(x, t) = 0. The BCs in time 
involve kink- like particle density profiles qi(x) and q2{x), 
at t = and t = r, respectively, separated by distance 
X = ct. Once the Hamilton's equations are solved, one 
can calculate the action S along the whole trajectory 
q(x,t), p(x,t) and evaluate the probability V(c): 

/OO r-T 
dx / dt[p(x,t)d t q-w] 
-oo JO 

/OO PT 
dx{AF+ dt[P(x,t)d t Q- W}}, (10) 
-oo JO 

where AF is the increment of the generating function 
density ([6]) between t — and t, and N = Klp/h I. 

Our crucial assumption, see also Ref. [18j| . is that 
the probability to observe, during a sufficiently long 
time r, a fluctuating front with a given front position 
X is determined by the action calculated along a TFS, 
Q = Q(x-ct) and P = P(x-ct) of Eqs. © and ©. This 
TFS solves the coupled ordinary differential equations 

Q" + cQ' + Q -Q 2 + 2QP-2Q 2 P = 0, (11) 
P" -cP' + P + P 2 - 2QP - 2QP 2 = . (12) 

subject to the proper BCs. Equations (fTTj) and (ITS]) pos- 
sess a conservation law: %q[Q{£,), P{Q] + Q'P' = const, 
where the constant is zero because of the BCs. 

For a TFS Eq. (|10|) can be simplified. As we will see, 
AF = for c > — 2. In this case we obtain 

-lnV(c)~KS = NTS, (13) 

with rescaled action accumulation rate 

/oo poo 
d£(cPQ' + W)= / d£,{Q-Q 2 )P 2 . (14) 
-oo J —oo 

We will also see that, for c < —2, it is the term PdtQ — W 
in Eq. (|10p that vanishes, and the only contribution to S 
comes from AF. 

The model A 2A and random walk, as described by 
the WKB theory, has two remarkable symmetries: 

1. Let Qi(x,t) — u(x,t) and Pi(x,t) — v(x,t) be 
a solution to Eqs. © and © obeying BCs A. Then 
Qi — —v(—x, —t) and Pi = —u(—x, — t) is also a so- 
lution obeying these BCs. Notice, that Q and P inter- 
changed. In particular, for TFSs we have Q\ = u(x — ct) 
and P\ — v{x — ct), and Q2 — —v(—x + ct) and Pi{x,t) — 
—u(—x + ct). Now, if there is a unique TFS (up to a shift 
in £), we obtain 

P(0 = -Q(to-Z), (15) 

that is Q and P for the same TFS are related. 

2. Let Qi(x,t) and Pi{x, t) be a solution to Eqs. © 
and © obeying BCs A. Then Q 2 = Pi(x,-t) + 1 and 
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P% = Qx(x, —t) — 1 is also a solution obeying these BCs. 
Q and P again interchanged! This yields a relation be- 
tween fluctuating TFSs with velocity c and — c: 



P(-c)(0 = Q(c) (£ + const) - 1. 



(16) 



In variables q and p this is simply qr_ c \(£) — <?( c )(£ + 
const): reversibility stemming from the detailed balance 
property of the microscopic model. As a result, the ac- 
tions for c < and — c > arc simply related: 



1(c) 



l(-c) 



(17) 



providing one more example of the Onsager-Machlup re- 
lation [23j. 

c = 0: Exact solution. How rare is the situation 
when the front motion is arrested by noise during a suf- 
ficiently long time r? For c = Eq. p^|) reduces to 
P(x) = Q(x) - 1, and the solution of Eqs. © and © 
with BCs A is elementary: 



Q(x) = (l + e x Y 



P(x) = -(l + e-*)" 



(18) 



where we have arbitrarily fixed the front position. Here 
AF = 0, and we can evaluate the action from Eqs. (p~3|) 
and dHJ). We obtain s = 1/3, so V(c = 0) ~ 
exp(— Nt/3). In the variables q and p the solution (fT£ 



becomes g(x) = (1 



and p(x) = — ln(l 



c < — 2: Exact solution. Another important obser- 
vation is that Q(x,t) — 1 solves the Hamilton's equation 
(|8]) and obeys the BCs B. As Q — qe~ p , we realize that 
p = hiq is an invariant manifold of Eqs. (j4]) and ([5]) for 
q and p. Elimination of p reduces each of Eqs. (j4|) and 
J3]) to d t q = — (q — q 2 + d%q), which is a time-reversed 
Eq. HU). The legitimate TFSs, with q(£ = -oo) = 1 and 
g(£ = oo) = 0, are those with q > 0. Therefore, we arrive 
a t <Z(c)(£) = Qo,-c(£ + const), for any c € (— oo, —2]. 

Now we can calculate the action for these fronts, using 
the first line of Eq. (flQ]l . As p — Inq, w = 0. The 
integral over t can be evaluated by parts, giving AF = 
f[q(x, r)] — f[q(x, 0)], where f(q) = q\nq—q. We evaluate 
AF on the TFS with a given c. The integration over x 
yields — cr, and we again arrive at Eq. (|13|) . with s = — c. 
The same result follows from the second line of Eq. ([TO]) , 
where the only contribution comes from AF. 

The origin of the simplicity of the case of c < —2 is in 
the equilibrium properties of the reversible microscopic 
model. In the discrete version of the model, the equi- 
librium multivariate distribution of the number of par- 
ticles in a finite system with N max lattice sites, is the 

N max >i . _ K 

product Poisson distribution [12j: P(n) = f| K * e , — . 



Assuming m 3> 1 and using the Stirling's formula, 
one obtains lnP(n) ~ — NV{q(x)}, where V{q(x)} — 
J dx {q\nq - q + 1) = J Q dxF(q,l), L = N max h, and 
we have returned to the continuous descriptions in q and 
x. Now one can consider, for L —¥ oo, a TFS moving with 



velocity c < — 2 and evaluate V{c) arriving at Eq. (fl3|) 
with s — —c. Note that, as the function s(c) = — c is 
not convex, TFSs may not correspond to the most likely 
particle density profiles for c < — 2. 

Numerics. For an arbitrary — 2 < c < 2 one needs to 
solve Eqs. (fTTj) and (fP2|) numerically. We used a shooting 
algorithm described elsewhere [l8j]. In view of Eqs. (fl())) 
and p7[) it would suffice to find the numerical solutions 
for < c < 2. Still, we computed several cases of c < 0, 
and verified this symmetry. One example of numerically 
found Q and P is shown in Fig. Q] Figure [2ji shows the 
numerically evaluated s, see Eq. (fT4)) . at different c. 
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FIG. 1: Numerical Q(£) and P(£) profiles for c = 1.5. 
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FIG. 2: (a) Symbols: the numerically found action accumula- 
tion rate s, see Eq. (| 14[) . versus the front velocity c. Straight 
line: the asymptote s = — c which becomes exact at c < —2. 
(b) Ins versus — c for c > 1.5. Symbols: numerical 

results. Dashed line: asymptote (|22[) . 



One can notice from Fig. [T]that the P-profile is shifted 
"downstream" compared with the Q-profile. The shift is 
equal to zero for c = 0, see Eq. (fTS)) , and increases with 
c. As c approaches 2, the shift becomes very large. That 
is, Q(£) first goes down to zero, while P ~ remains 



4 



almost unchanged. Then P goes down to —1 while Q 
is already close to zero. This salient feature calls for a 
perturbation theory that we will now present. 

2 - c <c 1: Nonlinear perturbation theory. At 
2 — c -C 1 the main contribution to the action comes from 
the leading edge of the front, £ >■ 1, where Q <C 1. As 
\P\ ~ 1 there, this theory must be nonlinear. This is in 
contrast to the front propagation into a metastable state, 
where a linear theory holds for c close to its mean-field 
value {18]. As one can see from the numerical solutions, 
\P\ <C 1 in the left region, whereas Q <C 1 in the right 
region. An analytical theory is possible because of the 
presence of a joint region, where \P\ <C 1 and Q « 1 
simultaneously. 

In the left region we can drop terms with P in Eq. (| 1 1|) 
and approximate by the mean-field solution <3o,c(£)j 
where 2 - c < 1. Now we linearize Eq. ([HI): P" - cP' + 
P — 2Q 0yC (^)P — 0. The solution decaying at £ — > — oo 
is Pi e f t {Q = PoQb. c (£,) eC ^ i with a yet unknown constant 
Po <C 1. In the right region we can drop terms with Q 
in Eq. flU): P" - cP' + P + P 2 = 0. Because of the 
symmetry (1151) . the solution is 



Pright(0 — ~ Qo,c(£o — £) • 



(19) 



with a yet unknown > 1. To find po and £o> we match 
the asymptotes, Pi e ft{Q = PrightiQ, in the joint region. 
Here both £ and £o — £ are much greater than 1, and 
we can use the leading-edge asymptote Qo,c{£, 3> 1) 



A(2 



-1/2 



e £ sin (\/2 — c £), where A ~ 1, and we 



have fixed the front position so that the cosine term is 
absent. The matching yields 



£o = 7r/V2 _ ^ + 1 + O (V2~c) 



Po 



(20) 



Now P(£) is known, in terms of the overlapping asymp- 
totes P ie ft and P r ight, for all f. By symmetry, 



Q(0 



Qo,c(0 , 



Co - e » i, 
£»1. 



(21) 



Now we calculate s from (|14p . The main contribution to 
the integral comes from the leading edge, £ 3> 1, where 
Q 2 can be neglected. Using Eqs. (fl9|) and (f2"Tj) and keep- 



ing the leading and subleading terms in i/2 — c in the 
exponent, we obtain after integration by parts 



3e 



Ql , 2 (C)dC- 0.0074 e" 



(22) 



where <3o,2(£) i s fixed by the demand that Qo,2(£ 3> 1) — 
A£e~£, whereas the e - ^ term is absent. 

In SUMMARY, we evaluated, for the microscopic model 
A «=2 2 A and random walk, the probability V(c) of ob- 
serving a fluctuating FKPP front moving with velocity 
c < 2: sec Eq. (JTSJ) and Fig. [2j For 2 - c < 1 we found 



This result holds only when there are many particles at 
the leading edge of the front, where £ ~ £o — 7r/\/2 — c. 
This demands c* — c ^> 27r 2 In -3 A, where c* = 2 — 
7r 2 In -2 iV. Somewhat surprisingly, c* coincides with the 
velocity of the mean-field FKPP front with account of 
discreteness of particles Q. In view of the strong in- 
equality c* — c ^> 27r 2 In -3 N, Eq. (|23|) may have a joint 
region of validity with the results of Brunet et al. [l5| . 
This issue demands a careful study [24|. We only note 
here that, by putting c = c* — 5c, where <5c <C In -2 A, the 
linear dependence on N in Eq. (|23l) cancels out, leading 
to a much slower logarithmic dependence, cf. Ref. 



15] 



Finally, we have found that Eq. (|23l) holds, up to a 
c-independent coefficient, for all pulled reaction fronts, 
where A — > 2A is the only first-order birth process. It 
also holds, close to the transition, for all sets of reactions 
belonging to the Directed Percolation universality class. 
The properly rescaled on-site Hamiltonian for this class 
of models is H (Q,P) = QP(P-Q + 1) 
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